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Understanding how institutional changes within academia may affect the overall potential of sci- 
ence requires a better quantitative representation of how careers evolve over time. Since knowledge 
spillovers, cumulative advantage, competition, and collaboration are distinctive features of the aca- 
demic profession, both the employment relationship and the procedures for assigning recognition 
and allocating funding should be designed to account for these factors. We study the annual pro- 
duction ni(t) of a given scientist i by analyzing longitudinal career data for 200 leading scientists 
and 100 assistant professors from the physics community. We compare our results with 21,156 
sports careers. Our empirical analysis of individual productivity dynamics shows that (i) there are 
increasing returns for the top individuals within the competitive cohort, and that (ii) the distribu- 
tion of production growth is a leptokurtic "tent-shaped" distribution that is remarkably symmetric. 
Our methodology is general, and we speculate that similar features appear in other disciplines 
where academic publication is essential and collaboration is a key feature. We introduce a model 
of proportional growth which reproduces these two observations, and additionally accounts for the 
significantly right-skewed distributions of career longevity and achievement in science. Using this 
theoretical model, we show that short-term contracts can amplify the effects of competition and 
uncertainty making careers more vulnerable to early termination, not necessarily due to lack of 
individual talent and persistence, but because of random negative production shocks. We show that 
fluctuations in scientific production are quantitatively related to a scientist's collaboration radius 
and team efficiency. 



Institutional change could alter the relationship be- 
tween Science and scientists as well as the longstanding 
patronage system in academia [TJ [5] . Some recent shifts 
in academia include the changing business structure of re- 
search universities [3], shifts in the labor supply-demand 
balance a bottleneck in the number of tenure track 
positions [5], and a related policy shift away from long- 
term contracts E] • Along these lines, significant fac- 
tors for consideration are the increasing range in research 
team size [7] , the economic organization required to fund 
and review collaborative research projects, and the evolv- 
ing definition of the role of the academic research profes- 
sor [3]. 

The role of individual performance metrics in career 
appraisal, in domains as diverse as sports [5J [5] , finance 
[TtJl ITT] . and academia, is increasing in this data rich 
age. In the case of academia, as the typical size of scien- 
tific collaborations increases [7] , the allocation of funding 
and the association of recognition at the varying scales 
of science (individual group ^=> institution [T^]) has 
become more complex. Indeed, scientific achievement is 
becoming increasingly linked to online visibility in a con- 
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siderable reputation tournament |13j . 

Here we seek to identify (i) quantitative patterns in the 
scientific career trajectory towards a better understand- 
ing of career dynamics and achievement [T4] - (20] , and (ii) 
how scientific production responds to policies concerning 
contract length. Using rich productivity data available 
at the level of single individuals, we analyze longitudinal 
career data keeping in mind the the roles of spillovers, 
group size, and career sustainability. Although our em- 
pirical analysis is limited to careers in physics, our ap- 
proach is general. We speculate that similar features de- 
scribe other disciplines where academic publication is a 
primary indicator and collaboration is a key feature. 

Specifically, we analyze production data for 300 physi- 
cists i = 1...300 who are distributed into 3 groups: (a) 
Group A corresponds to the 100 most cited physicists 
with average h-mdex (h) = 61 ± 21, (b) Group B cor- 
responds to 100 additional highly-cited physicists with 
(h) = 44 ± 15, and (c) Group C corresponds to 100 as- 
sistant professors in 50 U.S. physics departments with 
(h) = 15 ± 7. We define the annual production n.;(i) as 
the number of papers published by scientist i in year t 
of his/her career. We focus on academic careers from 
the physics community to approximately control for sig- 
nificant cross-disciplinary production variations. Using 
the same set of scientists, a companion study has ana- 
lyzed the rank-ordered citation distribution of each sci- 
entist with a focus on the statistical regularities underly- 
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ing publication impact [T7]. We provide further descrip- 
tion of the data and present a parallel analysis of 21,156 
sports careers in the Supporting Information Appendix 
(SI) text. 

We begin this paper with empirical analysis of longi- 
tudinal career data. Our empirical evidences serve as 
statistical benchmarks used in the final section where we 
develop a stochastic proportional growth model. In par- 
ticular, our model shows that a short-term appraisal sys- 
tem can result in a significant number of "sudden" early 
deaths due to unavoidable negative production shocks. 
This result is consistent with a Matthew Effect model 
[16j and recent academic career survival analysis |21j . 
which demonstrate how young careers can be stymied 
by the difficulty in overcoming early achievement barri- 
ers. Altogether, our results indicate that short-term con- 
tracts may increase the strength of the "rich-get-richer" 
mechanism in science [22 , 23J and may hinder the upward 
mobility of young scientists. 

I. RESULTS 
A. Scientific production and the career trajectory 

The academic career depends on many factors, such as 
cumulative advantage [THJ HH H21 H3] , the "sacred spark," 
[24l [25] , and other complex aspects of knowledge trans- 
fer manifest in our techno-social world [26]. To exem- 
plify this complexity, a recent case study on the impact 
trajectories of Nobel prize winners shows that "scientific 
career shocks" marked by the publication of an individ- 
ual's "magnum opus" work(s) can trigger future recogni- 
tion and reward, resembling the cascading dynamics of 
earthquakes [2"?] . 

We model the career trajectory as a sequence of scien- 
tific outputs which arrive at the variable rate rii(t). Since 
the reputation of a scientist is typically a cumulative rep- 
resentation of his/her contributions, we consider the cu- 
mulative production iVj(t) = 53t'=i n i(t') as a proxy for 
career achievement. Fig. [TjA) shows the cumulative pro- 
duction Ni (t) of six notable careers which display a tem- 
poral scaling relation Ni(i) ps Ait ai where cti is a scaling 
exponent that quantifies the career trajectory dynam- 
ics. The average and standard deviation of the on values 
calculated for each dataset are (on) = 1.42 ± 0.29 [A], 
1.44 ± 0.26 [B], and 1.30 ± 0.31 [C]. We justify this 2- 
parameter model in the SI Appendix text using scaling 
methods and data collapse. 

There are also numerous cases of iVj(t) which do not 
exhibit such regularity (see Fig. SI), but instead display 
marked non-stationarity and non-linearity arising from 
significant exogenous career shocks. Positive shocks, pos- 
sibly corresponding to just a single discovery, can spur 
significant productivity and reputation growth [24] |27] . 
Negative shocks, such as in the case of scientific fraud, 
can end the career rather suddenly. We also acknowledge 
that the end of the career is a difficult phase to analyze, 
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FIG. 1: Persistent accelerating career growth. (A) The ca- 
reer trajectory Ni(t) ~ t ai of six stellar careers from varying 
age cohorts. The cti value characterizes the career persis- 
tence, where careers with a > 1 are accelerating, an val- 
ues calculated using OLS regression in alphabetical order 
are: a = 1.25 ± 0.02, a = 1.72 ± 0.02, a = 1.62 ± 0.04, 
a = 1.23±0.02, a = 1.34±0.05, a = 1.35±0.04. (B) Denned 
in Eq. [lj, the average career trajectory (N'(t)) calculated 
from 100 individual Ni(t) in each dataset demonstrates ro- 
bust accelerating career growth within each cohort. We use 
the normalized career trajectory Nl(t) in order to aggregate 
Ni(t) with varying publication rates (rii). As a result, the ag- 
gregate scaling exponent a quantifies the acceleration of the 
typical career over time, independent of (rii). For the scien- 
tific careers, we calculate a values: 1.28±0.01 [A], 1.31±0.01 
[B], and 1.15 ± 0.02 [C]. These values are all significantly 
greater than unity, a > 1, indicating that cumulative advan- 
tage in science is closely related to knowledge and production 
spillovers. We calculate a using OLS regression and plot the 
corresponding best-fit lines (dashed) for each dataset. 



since such an event can occur quite abruptly, and so our 
analysis is mainly concerned with the growth phase and 
not the termination phase. 

In order to analyze the average properties of Ni(t) for 
all 300 scientists in our sample, we define the normalized 
trajectory Nl(t) = Ni(t)/(rii). The quantity (rii) is the 
average annual production of author i, with N[(Li) = L t 
by construction (Li corresponds to the career length of 
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individual i). Fig. [T[B) shows the characteristic produc- 
tion trajectory obtained by averaging together the 100 
N-(t) belonging to each dataset, 

^ v (t)) - \ (n t ) I ~ ioo ^ (m) ■ [L) 

The standard deviation a(N'(t)) shown in Fig. S2(B) 
begins to decrease after roughly 20 years for dataset [A] 
and [B] scientists. Over this horizon, the stochastic ar- 
rival of career shocks can significantly alter the career 
trajectory [HI Hfl W\ [25] • Each N{(t) exhibits robust 
scaling corresponding to the scaling law (N'(t)} ~ t a . 
This regularity reflects the abundance of of careers with 
oti > 1 corresponding to accelerated career growth. This 
acceleration is consistent with increasing returns arising 
from knowledge and production spillovers. 



B. Fluctuations in scientific output over the 
academic career 

Individuals are constantly entering and exiting the 
professional market, with birth and death rates depend- 
ing on complex economic and institutional factors. Due 
to competition, decisions and performance at the early 
stages of the career can have long lasting consequences 
[TBI |2T)] . To better understand career uncertainty por- 
trayed by the common saying "publish or perish" [30] . 
we analyze the outcome fluctuation 

n(t) = m(t) - m(t - At) (2) 

of career i in year t over the time interval At — 1 year. 
Fig. [2jA) and (B) show the unconditional pdf of r values 
which arc leptokurtic but remarkably symmetric, illus- 
trating the endogenous frequencies of positive and nega- 
tive output growth. Output fluctuations arise naturally 
from the lulls and bursts in both the mental and physical 
capabilities of humans [3TJ|32]. Moreover, the statistical 
regularities in the annual production change distribution 
indicate a striking resemblance to the growth rate distri- 
bution of countries, firms, and universities [3"3l l3~i] . 

To better account for individual growth factors, we 
next define the normalized production change 

ri(t) = [ri(t)-{ri)]/<Ti(r) (3) 

which is measured in units of the fluctuation scale <7j(r) 
unique to each career. We measure the average (j-j) and 
the standard deviation <Xj(r) of each career using the first 
Li available years for each scientist i. r'^t) is a better 
measure for comparing career uncertainty, since individ- 
uals have production factors that depend on the type 
of research, the size of the collaboration team, and the 
position within the team. Fig. [2jC) show that P(r'), 
the probability density function (pdf) of r' measured in 
units of standard deviation, is well approximated by a 




FIG. 2: Empirical evidence for the proportional growth 
model of career production. (A) Probability density func- 
tion (pdf) of the annual production change r in the number 
of papers published over a At = 1 year period. In the bulk of 
each P(r), the growth distribution is approximately double- 
exponential (Laplace). (B) To test the stability of the distri- 
bution over career trajectory subintervals, we separate r»(t) 
values into 5 non-overlapping 10-year periods and verify the 
stability of the Laplace P(r). For each P(r), we also plot the 
corresponding Laplace distribution (solid line) with standard 
deviation a and mean /i « calculated using the maximum 
likelihood estimator method. To improve graphical clarity, 
we vertically offset each P(r) by a constant factor. For vi- 
sual comparison, we also plot a Normal distribution (dashed 
black curve) with a = 1 which instead decays parabolically 
on the log-linear axes. (C) Accounting for individual produc- 
tion factors by using the normalized production change r', 
the resulting pdfs P(r') collapse onto a Gaussian distribution 
with unit variance. Deviations in the tails likely correspond 
to extreme "career shocks." (D) The cumulative distribution 
CDF(X > Si) is exponential, indicating that the uncondi- 
tional distributions P(r) in (A) and (B) follow from an expo- 
nential mixing of conditional Gaussian distributions P(r\Si). 



Gaussian distribution with unit variance. The data col- 
lapse of each P(r') onto the predicted Gaussian distribu- 
tion (solid green curve) indicates that individual output 
fluctuations are consistent with a proportional growth 
model. We note that the remaining deviations in the 
tails for \r'\ > 3 are likely signatures of the exogenous 
career shocks that are not accounted for by an endoge- 
nous proportional growth model. 

The ability to collaborate on large projects, both in 
close working teams and in extreme examples as remote 
agents (i.e. Wikipedia [35]), is one of the foremost prop- 
erties of human society. In science, the ability to attract 
future opportunities is strongly related to production and 
knowledge spillovers [35H3BJ[37] that are facilitated by the 
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collaboration network [71 [T^l I58H4"2] . Indeed, there is a 
tipping point in a scientific career that occurs when a sci- 
entist's knowledge investment reaches a critical mass that 
can sustain production over a long horizon, and when 
a scientist becomes an attractor (as opposed to a pur- 
suer) of new collaboration/production opportunities. To 
account for collaboration, we calculate for each author 
the number fc, (t) of distinct coauthors per year and then 
define his/her collaboration radius Si as the median of 
the set of his/her fcj(i) values, Si = Med[ki(t)]. We use 
the median instead of the average (ki(t)) since extremely 
large ki (t) values can occur in specific fields such as high- 
energy physics and astronomy. 

Given the complex scientific coauthorship network, 
we ask the question: what is the typical number of 
unique coauthors per year? Fig. (2TD) shows the cu- 
mulative distribution function CDF[Si) of Si values for 
each data set. The approximately linear form on log- 
linear axes indicates that Si is exponentially distributed, 
P(Si) ~ exp[-A5i]. We calculate A = 0.15 ± 0.01 [A], 
A = 0.11 ± 0.01 [B], and A = 0.11 ± 0.01 [C]. The ex- 
ponential size distribution has been shown to emerge 
in complex systems where linear preferential attachment 
governs the acquisition of new opportunities |43j . This 
result shows that the leptokurtic "tent-shaped" distribu- 
tion P(r) in Fig. [2] follows from the exponential mixing 
of heterogenous conditional Gaussian distributions [H] . 

The exponential mixture of Gaussians decomposes the 
unconditional distribution P{r) into a mixture of condi- 
tional Gaussian distributions 



P{r\Si) = exp{~r 2 /2VSt}/^2nVSt , 



(4) 



each with a fluctuation scale &i(r) depending on Si by 
the scaling relation 

af(r) « VSf . (5) 
Hence, the mixture is parameterized by ip 

P^{r)= P{r\S)P{S)dS^Y. P ^ r \ S ^ P ^) ■ ( 6 ) 
J o i=l 

The independent case ip = results in a Gaussian P^(r) 
and the linear case ip = 1 results in a Laplace (double- 
exponential) P^(r). See the SI Appendix text and ref. 
[4"4"j for further discussion of the ip dependence of P$ (r) . 



C. The size- variance relation and group efficiency 

The values of ip for scientific and athletic careers follow 
from the different combination of physical and intellec- 
tual inputs that enter the production function for the 
two distinct professions. Academic knowledge is typi- 
cally a non-rival good, and so knowledge-intensive pro- 
fessions are characterized by spillovers, both over time 
and across collaborations [351 13Z] , consistent with a. t > 1 
and ip > 0. Interestingly, Azoulay et al. show evidence 
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FIG. 3: Quantitative relations between career growth, career 
risk, and collaboration efficiency. The fluctuations in produc- 
tion reflect the unpredictable horizon of "career shocks" which 
can affect the ability of a scientists to access new creative op- 
portunities. (A) Relation between average annual production 
(rn) and collaboration radius Si = M ed[ki] shows a decreasing 
marginal output per collaborator as demonstrated by sublin- 
ear ip < 1. Interestingly, dataset [A] scientists have on average 
a larger output-to-input efficiency. (B) The production fluc- 
tuation scale (Ji{r) is a quantitative measure for uncertainty 
in academic careers, with scaling relation Oi(r) ~ Sf^ 2 . (C) 
Management, coordination, and training inefficiencies can re- 
sult in a 7 < 1 corresponding to a decreasing marginal re- 
turn with each additional coauthor input. The significantly 
larger 7 value for dataset [A] scientists seems to suggest that 
managerial abilities related to output efficiency is a common 
attribute of top scientists. 
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for production spillovers in the 5-8% decrease in output 
by scientists who were close collaborators with a "super- 
star" scientists who died suddenly [25] . 

We now formalize the quantitative link between scien- 
tific collaboration 38, 39] and career growth given by the 
size-variance scaling relation in Eq. [5] visualized in the 
scatter plot in Fig. [3](B). Using ordinary least squares 
(OLS) regression of the data on log-log scale, we cal- 
culate ip/2 ps 0.40 ± 0.03 (R = 0.77) for datasct [A], 
ip/2 ps 0.22 ±0.04 (R = 0.51) [B], and tjj/2 ps 0.26 ±0.05 
(R = 0.45) [C]. Interdependent tasks characteristic of 
group collaborations typically involve partially overlap- 
ping efforts. Hence, the empirical ip values are signifi- 
cantly less than the value ip = 1 that one would expect 
from the sum of Si independent random variables with 
approximately equal variance V. Collectively, these em- 
pirical evidences serve as coherent motivations for the the 
preferential capture growth model that we propose in the 
following section. 

Alternatively, it is also possible to estimate tp using 
the relation between the average annual production (rij) 
and the collaboration radius Sj. The input-output re- 
lation (n,-) ~ Sf quantifies the collaboration efficiency, 
with -0 = 0.74 ± 0.04 (R = 0.87) for dataset [A] and 
tp = 0.25±0.04 (R = 0.37) for dataset [B]. If the autocor- 
relation between sequential production values rii(t) and 
rii(t + 1) is relatively small, then we expect the scaling 
exponents calculated for (n,) and of(r) to be approxi- 
mately equal. This result follows from considering Ti(t) 
as the convolution of an underlying production distribu- 
tion Pj(n) for each scientist that is approximately stable. 
Interestingly, the larger -0 values calculated for dataset 
[A] scientists suggests that prestige is related to the in- 
creasing returns in the scientific production function |45j . 

Next we use an alternative method to estimate the 
annual collaboration efficiency by relating the number 
of publications rii(t) in a given year to the number of 
distinct coauthors fcj(t) over the same year. We use a 
single-factor production function, 



rii(t) ps qi[ki(t)} 



(7) 



to quantify the relation between output and labor in- 
puts with a scaling exponent 7.;. We estimate qi and 
7i for each author using OLS regression, and define the 
normalized output measure Qi oc n.j(i) / ki(t) 7i using the 
best-fit qi and 7* values calculated for each scientist i. 
Fig. [3jC) shows the efficiency parameter 7 calculated 
by aggregating all careers in each dataset, and indicates 
that this aggregate 7 is approximately equal to the av- 
erage (tj) calculated from the 7, values in each career 
dataset: 7 = 0.68 ± 0.01 [A], 7 = 0.52 ± 0.01 [B], and 
7 = 0.51 ± 0.02 [C]. Furthermore, the ip and 7 values are 
approximately equal, which is not surprising, since both 
scaling exponents are efficiency measures that relate the 
scaling relation of output rii(t) per input ki(t). 



D. A Proportional growth model for scientific 
output 

We develop a stochastic model as a heuristic tool to 
better understand the effects of long-term versus short- 
term contracts. In this competition model, opportunities 
(i.e. new scientific publications) are captured according 
to a general mechanism whereby the capture rate Vi(t) 
depends on the appraisal Wi(t) of an individual's record 
of achievement over a prescribed history. We define the 
appraisal to be an exponentially weighted average over a 
given individual's history of production 



t-i 



Oi(t) = J2 ih(t-At)e 



-cAt 



(8) 



At=l 



which is characterized by the appraisal horizon 1/c. We 
use the value c = to represent a long-term appraisal 
(tenure) system and a value c ^> 1 to represent a short- 
term appraisal system. Each agent i = simultane- 
ously attracts new opportunities at a rate 



wi{ty 



(9) 



until all P opportunities for a given period t are cap- 
tured. We assume that each agent has the production 
potential of one unit per period, and so the total number 
of opportunities distributed per period P is equal to the 
number of competing agents, P = I. 

We use Monte Carlo (MC) simulation to analyze this 
2-parameter model over the course of t = 1...T sequential 
periods. In each production period (i.e. representing a 
characteristic time to publication), a fixed number of P 
production units are captured by the competing agents. 
At the end of each period, we update each Wi (t) and then 
proceed to simulate the next preferential capture period 
t+1. Since Vi(t) depends on the relative achievements of 
every agent, the relative competitive advantage of one in- 
dividual over another is determined by the parameter ir. 
In the SI Appendix text we elaborate in more detail the 
results of our simulation of synthetic careers dynamics. 
We vary ir and c for a labor force of size I = 1000 and 
maximum lifetime T = 100 periods as a representative 
size and duration of a real labor cohort. Our results are 
general, and for sufficiently large system size, the quali- 
tative features of the results do not depend significantly 
on the choice of I or T. 

The case with n = corresponds to a random capture 
model that has (i) no appraisal and (ii) no preferential 
capture. Hence, in this null model, opportunities are 
captured at a Poisson rate X p = 1 per period. The results 
of this model (see Fig. S13) shows that almost all careers 
obtain the maximum career length T with a typical career 
trajectory exponent {on) ps 1. Comparing to simulations 
with 7r > and c > 0, the null model is similar to a 
"long-term" appraisal system (c — > 0) with sub-linear 
preferential capture (n < 1). In such systems, the long- 
term appraisal timescale averages out fluctuations, and 
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FIG. 4: Monte Carlo simulation of the linear preferential capture model (tt = 1) for varying contract length parametrized by c. 
We plot the probability distributions for (i) N,, the total number of opportunities captured by the end period T, (ii) the growth 
acceleration exponent a%, (iii) the single period growth fluctuation ri(t) including for comparison the Laplace (solid green) 
and Gaussian (dashed red) best-fit distributions calculated using the respective MLE estimator, and (iv) the career longevity 
Li defined as the time difference between an agent's first and last captured opportunity. Results for c — > systems shows 
that for a "long-term appraisal" scenario careers are less vulnerable to low-production phases, and as a result, most agents 
sustain production throughout the career. Conversely, results for c > 1 systems show that for a "short-term appraisal" scenario 
the labor system is driven by fluctuations that can cause career "sudden death" for a large fraction of the population. In this 
short-term appraisal model, there are typically a small number of agents who are able to capture the majority of the production 
opportunities with remarkably accelerating career growth reflected by significantly large oti > 1. Thus, a few "lucky" agents are 
able to survive the initial fluctuations and end up dominating the system. In the SI text and Figs. S12-S16, we further show 
that systems with increased levels of competition (n > I) mimic systems with short term contracts, resulting in productivity 
"death traps" whereby most careers stagnate and terminate early. 



so careers are significantly less vulnerable to periods of 
low production and hence more sustainable since they are 
not determined primarily by early career fluctuations. 

However, as 7r increases, the strength of competitive 
advantage in the system increases, and so some careers 
are "squeezed out" by the larger more dominant careers. 
This effect is compounded by short-term appraisal cor- 
responding to c « 1. In such systems with super- linear 
capture rates and/or relatively large c, most individuals 
experience "sudden death" termination relatively early in 
the career. Meanwhile, a small number of "stars" survive 
the initial selection process, which is governed primarily 
by random chance, and dominate the system. 

We found drastically different lifetime distributions 
when we varied the appraisal (contract) length (see Figs. 
SI2 - S16). In the case of linear preferential capture 
with a long-term appraisal system c = 0, we find that 
10% of the labor population terminates before reaching 



career age 0.94T (where T is the maximum career length 
or "retirement age" ) , and only 25% of the labor popula- 
tion terminates before reaching career age 0.98T. On the 
contrary, in a short-term appraisal system with c = 1, 
we find that 10% of the labor population terminates be- 
fore reaching age 0.01T, and 25% of the labor population 
dies before reaching age 0.02T (see Table SI). Hence, in 
model short contract systems, the longevity, output, and 
impact of careers are largely determined by fluctuations 
and not by persistence. 

Fig. [3] shows the MC results for n = 1. For c > 1 we 
observe a drastic shift in the career longevity distribution 
P(L), which becomes heavily right-skewed with most ca- 
reers terminating extremely early This observation is 
consistent with the predictions of an analytically solv- 
able Matthew effect model [16] which demonstrates that 
many careers have difficulty making forward progress due 
to the relative disadvantage associated with early career 
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inexperience. However, due to the nature of zero-sum 
competition, there are a few "big winners" who survive 
for the entire duration T and who acquire a majority of 
the opportunities allocated during the evolution of the 
system. Quantitatively, the distribution P(N) becomes 
extremely heavy-tailed due to agents with a > 2 corre- 
sponding to extreme accelerating career growth. Despite 
the fact that all the agents are endowed initially with 
the same production potential, some agents emerge as 
superstars following stochastic fluctuations at relatively 
early stages of the career, thus reaping the full benefits 
of cumulative advantage. 



II. DISCUSSION 

An ongoing debate involving academics, university ad- 
ministration, and educational policy makers concerns 
the definition of professorship and the case for lifetime 
tenure, as changes in the economics of university growth 
have now placed tenure under the review process El E] . 
Critics of tenure argue that tenure places too much fi- 
nancial risk burden on the modern competitive research 
university and diminishes the ability to adapt to shift- 
ing economic, employment, and scientific markets. To 
address these changes, universities and other research in- 
stitutes have shifted away from tenure at all levels of 
academia in the last thirty years towards meeting staff 
needs with short-term and non-tenure track positions E] . 

For knowledge intensive domains, production is char- 
acterized by long-term spillovers both through time and 
through the knowledge network of associated ideas and 
agents. A potential drawback of professions designed 
around short-term contracts is that there is an implicit 
expectation of sustained annual production that effec- 
tively discounts the cumulative achievements of the in- 
dividual. Consequently, there is a possibility that short- 
term contracts may reduce the incentives for a young sci- 
entist to invest in human and social capital accumulation. 
Moreover, we highlight the importance of an employment 
relationship that is able to combine positive competitive 
pressure with adequate safeguards to protect against ca- 
reer hazards and endogenous production uncertainty an 
individual is likely to encounter in his/her career. 

In an attempt to render a more objective review pro- 



cess for tenure and other lifetime achievement awards, 
quantitative measures for scientific publication impact 
are increasing in use and variety [T?H2T)l EMI 1271 l4"rJI |4"7) . 
However, many quantifiable benchmarks such as the h- 
index [T7| do not take into account collaboration size 
or discipline specific factors. Measures for the compari- 
son of scientific achievement should at least account for 
variable collaboration, publication, and citation factors 
[Hfl US E7] . Hence, such open problems call for further 
research into the quantitative aspects of scientific output 
using comprehensive longitudinal data for not just the 
extremely prolific scientists, but the entire labor force. 

Current scientific trends indicate that there will be 
further increases in typical team sizes that will forward 
the emergent complexity arising from group dynamics 
[T^l • There is an increasing need for individ- 
ual/group production measures, such as the output mea- 
sure Q, following from Eq. [7 , which accounts for 
group efficiency factors. Normalized production mea- 
sures which account for coauthorship factors have been 
proposed in [191 146] , but the measures proposed therein 
do not account for the variations in team productivity. 

The complexity of large collaborations raises open 
questions concerning scientific productivity and the or- 
ganization of teams. We measure a decreasing marginal 
returns 7 < 1 with increasing group size which identi- 
fies the importance of team management. A theory of 
labor productivity can help improve our understanding 
of institutional growth, for organizations ranging in size 
from scientific collaborations to universities, firms, and 
countries E3 El E 
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I. DATA 

To test the intriguing possibility that competition leads to common growth patterns in complex systems of arbitrary 
size S, we analyze the production dynamics of two professions that are dissimilar in many regards, but share the 
common underlying driving force of competition for limited resources. In order to establish empirical facts that we 
believe are independent of the details of a given competitive profession, we analyze a large dataset of production rii(t) 
values and corresponding growth fluctuation rj(t) = n^t) — i%i(t — 1) values. We define the appropriate measures 
for rii(t) to be (a) the annual number of papers published by scientist i and (b) the seasonal performance metrics 
of professional athlete i. While these two professions both display a high level of competition, they differ in their 
employment term structure and salary scale. In the case of academia, the tenure system rewards high performance 
levels with lifelong employment (tenure). In contrast, professional sports are characterized by relatively short contracts 
that emphasize continued performance over a shorter time frame and thereby exploit the high levels of athletic prowess 
in a player's peak years. The large number of careers in these two professions readily lend themselves to quantitative 
analysis because the data that quantify the career production trajectory are precisely defined and comprehensive 
throughout an individual's entire career. Furthermore, because of the generic nature of competition, we use these two 
distinct professions to compare and contrast the distribution of career impact measures across a cohort of competitors. 
The datasets we analyze are: 

I : Academia: 

We analyze the publication careers of 300 physicists which we categorize in 3 subsets each consisting of 100 
individuals: 

(A) Dataset A corresponds to the 100 most-cited physicists according to the citation shares metric [TH] 
(with average /i-index (h) — 61 ± 21). These 100 careers constitute 3,951 ri{t) values. 

(B) Dataset B corresponds to the 100 other "control" scientists, taken approximately randomly from the 
same physics database (with average /i-index (h) = 44 ± 15). In the selection process for dataset B, we 
only consider scientists who have published between 10 and 50 articles in PRL over the 50-year period 
1958-2008. These 100 careers constitute 3,534 r^t) values. 

(C) Dataset C corresponds to 100 Assistant Professors (with average h-'mdex (h) = 15 ±7), where we select 
two physicists from each of the top-50 U.S Physics & Astronomy Departments (according to the U.S. 
News rankings). These Asst. Profs, are assumed to be early in their career and relatively accomplished 
given the difficulty in obtaining such a position in any given university. These 100 careers constitute 
1,050 Ti{t) values. 
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In order to control for discipline-specific citation patterns, we select individuals in dataset A and B from set of all 
scientists who have published in Physical Review Letters (PRL) over the 50-year period 1958-2008. As a measure of 
output, we define rii(t) as the number of papers published in year t of the career of individual i, where year t = 1 
corresponds to the year of the first publication on record for author i. We downloaded the complete publication 
records of the scientists in datasets A and B from ISI Web of Science ( jhttp : //www, i siknowledge . com/) in Jan. 
2010, and we downloaded the complete publication records of the scientists in dataset C from ISI Web of Science in 
Oct. 2010. We used the "Distinct Author Sets" function provided by ISI in order to increase the likelihood that only 
papers published by each given author are analyzed. 

II : Major League Baseball (MLB): 

We analyze 17,292 baseball players over the 90-year period 1920-2009 using comprehensive league data obtained 
from Sean Lahman's Baseball Archive accessed at http://baseballl.com/index.php. We separate the career 
data into two distinct subsets: non-pitchers (players not on record as having pitched during a game) and pitchers. 

(A) For non-pitchers, we analyze two batting metrics: an "opportunity metric" - at-bats (AB), and a 
"success" metric - hits (H). Together, these 8,993 careers constitute 43,043 ri(t) values. 

(B) For pitchers, we analyze two pitching metrics: an "opportunity metric" - innings-pitched measured in 
outs (IPO), and a "success" metric - strikeouts (K). Together, these 8,299 careers constitute 33,965 
Ti{t) values. 

Ill : National Basketball Association (NBA): 

We analyze 3,864 basketball careers, constituting 15,316 r,(f) values, over the 63-year period 
1946-2008 using data obtained from Data Base Sports Basketball Archive accessed at 
|nttp : //www . d atabasebasketball . com/ . We analyze two player metrics: 

(A) an "opportunity metric" - minutes played (Min.), and 

(B) a "success" metric - points scored (Pts.) 

Since sports careers typically peak for athletes around age 30, we account for a time-dependent career trajectory 
which is dominant in most sports careers by "detrended" the measures for career growth fluctuations. In the case 
where we do not account for a individual fluctuation scale, 

Ri= [n(t)-r{t)]/a(t) . (SI) 

In this case we detrend with respect to the average production difference f(t) and the standard deviation of production 
difference o~(t) which are calculated using all careers from a given sports league, conditional on the career year t. 

In the case where we do account for individual variations, we first define Zi(t) = (ri(t) — (r"i})/<7j to be normalized 
with respect to the individual career scales (ri) and Oi which are the average and standard deviation of the production 
change of athlete career i. Then we define the detrended growth rate as 

^ = [Zi(t) - {z{t))]/a z{t) , (S2) 

where in this case we detrend with respect to the average (z(t)) and standard deviation <J z (t) calculated by collecting 
all Zi(t) values for a given career year t. This detrending better accounts for the relatively strong time-dependent 
growth patterns in sports. 

In this section we analyze the annual production of scientists measured as the number of papers published rjj(i) 
over the period of a year. Using this measure does not account for the variability in the length of production, say in 
the number of pages, nor does it account for the impact of the paper, a quantity commonly approximated by a paper's 
citation number. Instead, we consider a simple definition that a scientific product is a final output of a collection of 
inputs. Furthermore, in science it is assumed that the peer review process establishes a quality threshold so that only 
manuscripts above a certain quality and novelty standard can be published and incorporated into the scientific body 
of knowledge. 

Prior theories of scientific production have also used the number of publications as a proxy for scientific output. 
In particular, the Shockely model |14j proposed a simple multiplicative factor model for the production rij(t) which 
predicts a log-normal distribution for P{n). An alternative null model for rii(t) is the Poisson process, which assumes 
that each individual is endowed with a rate parameter uj related to an individual's production factors. This model 
predicts a Poisson distribution for P{n). However, a shortfall of these models is that multiplicative parameters in 
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the Shockley model and the rate parameter ui are difficult to measure, especially if the set of individuals span a large 
range of production factors, and moreover, if the careers are non-stationary. 



Fig. S8 shows the unconditional probability distribution P(n) calculated by aggregating all Jij(t) values for all 
scientists and all years into an aggregate dataset. Naively, the distributions are well- fit by the Log-normal distribution, 
and so there is an apparent agreement with the multiplicative factor Shockley model. However, the distribution 
P(n) = E t =iP(n\Si) is the aggregate distribution constructed from 100 individual career trajectories n z (t), each 
with varying size Si- Indeed, we demonstrate in Figs, fl] and [Sl| to be non-linear, with time-dependent residuals 
around the moving average. Hence, it is not possible from the unconditional pdf P(n) to determine if the process 
underlying scientific production corresponds to a simple multiplicative process or a Poisson process. 

In order to better account for the variable size Si of each career which affects the rate at which an individual is 
able to capture publication opportunities, we plot in Fig. |S7| the pdf of the normalized output 

* - if ■ *»> 

We calculate the normalization factor fi(k) = qi[ki{t)] li for each individual i by estimating the parameters qi and 7; 
for each scientist i from the single-factor model 

n % = qi kf . (S4) 

where n,(t) is the annual production in year t and ki(t) is the total number of distinct coauthors in year t. Hence, Qi 
represents the production factor above Q > 1 or below Q < 1 what would be expected from the author i given the fact 
that he/she had additional inputs from ki(t) — 1 individuals that year. This model assumes that the major component 
contributing to production is the collaboration degree k of the research output, and also assumes that the input of 
each coauthor contributes equally to the final output. Clearly, these assumptions neglect some important idiosyncratic 
details affecting scientific publication, but given the incomplete information associated with every publication, it is 
a decent approximation. We estimate q^ and 7; by performing a linear regression of log rii and log fcj using the first 
Li years of each career, neglecting years with rij = 0. We use Li = 35 years for dataset [A] and [B] scientists, and 
Li = 10 years for dataset [C] scientists. 

In Fig. [3jc) we approximate 7 using all n(t) within each dataset with k < 50, and performing a regression of the 
model 

In n = In q + 7m k + e (S5) 

to estimate 7, where e is the residual due to other unaccounted production factors. For each dataset we find that the 
aggregate efficiency parameter 7 is approximately equal to the average (7^) calculated from the 100 7, values in each 
career dataset: 7 = 0.68 ± 0.01 [A], 7 = 0.52 ± 0.01 [B], and 7 = 0.51 ± 0.02 [C]. Furthermore, the ^7 since the 
size-variance scaling parameter ip is also an efficiency measure that relates the scaling of output n to input k. 

As a result of this analysis, we quantify the scaling exponent 7 < 1 of the decreasing marginal returns in the scientific 
production function for projects with k < 50. This likely stems from the inefficient management costs associated with 
large group collaborations which typically manifest in a larger production timescale. In fact, for years with k > 50 
coauthors, scientific output shows decreasing returns to scale. Interestingly, the star scientists in dataset [A] display 
significantly larger efficiency, quantitatively showing the importance of management skills in scientific success. 

The normalized production values are normalized to units of "expected production" conditional on the ki inputs 
for author i. We aggregate all data from each dataset and show in Fig. [ST] that the Q values are well-described by 
the Gamma distribution 

TO) = Q--^f (SO) 

where m is the shape parameter and 9 is the scale parameter. Surprisingly, we find that dataset [A] and [B] have 
approximately equal Gamma parameters, indicating that besides their production efficiency, top scientists are virtually 
indistinguishable with average normalized output (Q) = m6 > 1. For each dataset we calculate the Gamma parameters 
using the maximum likelihood estimator method: m — 5.45 and 9 — 0.21 [A], m — 5.60 and 9 — 0.20 [B], and 
m = 7.00 and 9 — 0.15 [C]. We leave it as an open question to determine why the Gamma distribution describes so 
well the production statistics. We ponder the intriguing possibility that the stochastic dynamics underlying individual 
production corresponds to an increasing Levy process with variable jump length which is known to produce a Gamma 
distribution. 
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II. QUANTIFYING THE CAREER TRAJECTORY 



The reputation of an individual is typically cumulative, based on the total sum of achievements, which we approx- 
imate by the cumulative output Ni(t) (e.i. number of papers published by year t). In Figs, [ljand SI we plot -/Vj(i) 



for several individuals. The careers presented in Fig. [T]are more linear, indicating quantifiable career trajectory that 
has the approximate form 

t 

Ni(t) = n i(t') « A i tai . t < T i ( S7 ) 
{'=1 

where rii(t) are the number of papers in year t of the scientist's career which begins with t = 1 in the year of his/her 
first publication, and begins to decline around time Ti which is the time horizon over which the scaling regularity holds 
before termination and aging effects begin to dominate the career. In our analysis of academic career trajectories 
Ni(t), we only analyze Nj(t) for t < 40 years in order to account for such termination affects. 

The smooth career trajectories which appear as a linear curve when plotted on log-log scale are characterized by an 
amplitude parameter Ai and a scaling exponent a%. However, as indicated by Fig. |S1| there are also non-stationary 
Ni(t) which are dominated by "career shocks" that significantly alter the career trajectory. Such career shocks have 
been demonstrated using publication impact measures (e.i. citations, and h-index sequences) [3UJ [2H [S7J , and here 
we show that they even occur at the more fundamental level of individual production dynamics. 

In order to analyze the characteristic properties of iVj(i) for all 300 scientists analyzed, we define the normalized 
trajectory N-(t) = Nj(t)/{ni), where (n,i(t)) is the average annual production rate of author i, and so by construction 



Nl(Li) = Li. Fig. S2 A) shows the characteristic production trajectory obtained by averaging the 100 individual 



N-(t) for each dataset, 

ID!) 

100 



i=l 



The standard deviation a(N'(t)) is shown in Fig. |S2[ B), which has a broad peak that is a likely signature of 
career shocks that can significantly alter the career trajectory. The characteristic trajectory for each dataset are 
well-approximated by the scaling relation 

(N'(t)) ~ r (S9) 

with characteristic scaling exponents a > 1 that are significantly greater than unity: a = 1.28 ± 0.01 for Dataset A, 
a = 1.31 ±0.01 for Dataset B, and a = 1.15 ±0.02 for Dataset C. This fact implies that there is a significant cumulate 



advantage in scientific careers which allows for the career trajectory to be accelerating. In Fig. S2 'C) and S2 D) we 
plot the analogous (N'(t)) curves for professional sports metrics, where for this profession, a s» 1 for all measures 
analyzed. This quantitative feature is likely due to the fact that annual production in professional sports is capped 
by the limited number of opportunities provided by a season, whereas in academics, the number of publications a 
scientist can publish is in principle unlimited. Also, in more labour-intensive activities are likely to experience smaller 
returns since physical labor is non-cumulative with less spillover through time. 



In Fig. S3 we plot each individual career trajectory using the rescaled time t\ = t ai as an additional visual test of 



the scaling model given by Eq. S7 We show that on average, all curves i = 1..300 approximately collapse onto the 
expected curve Ni(t)/Ai = t', where the residual difference ej(t') = Ni(t)/Ai ~ t' are likely due to career shocks of 
various magnitudes. We plot the average and standard deviation of each set of 100 Ni(i)/Ai curves which show that 
most of the shocks Ci{t'), with some significant exceptions, lie within the la standard deviation denoted by the error 



bars. In Fig. S4 we plot the probability distributions P(aj) for each academic dataset. For each dataset, the average 



value (ai) is in good agreement with a, the scaling parameter calculated for the corresponding trajectory (N'(t)). 
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III. EXPONENTIAL MIXING OF GAUSSIANS 



The idea that entities are independent and identically distributed is an unrealistic assumption commonly made in 
analyses of complex systems. The unconditional pdf P(r) is commonly analyzed in empirical studies where insufficient 
data are present to define normalized r\ measures for each sample constituent i. Nevertheless, when modeling the 
evolution of complex based on empirical data corresponding to distinct subunits (such as individual careers, companies, 
or nation regions), unconditional quantities that account for variations in underlying production factors should be 
used. 

In the case of scientific output, there are many production factors that combine together and determine the amount 
of human efforts needed to produce a unit of production. In general, consider the value faj of individual i corresponding 
to his/her relative abilities in the production factor j = I... J corresponding to a variety of attributes: knowledge, 
genius, persistence, reputation, mental and physical health, communication skills, organization skills, and access to 
technology, equipment and data, etc. In this study, we compare scientists who publish in similar journals. Still, the 
scientific input required for each scientific output can vary by a large amount, largely depending on the technology 
needed to perform the analysis, ranging from particle accelerators to just a pencil and paper. 

In a very generalized representation, an unconditional distributions P(r), such as shown in Fig. [2]ja-d) for production 
change r, may follow from a mixture of conditional Gaussian distributions P(r\Si) 

P i> (r) = P(r\S)P(S)dS^y2p(r\S l )P(S l ) . (S10) 
J o <=1 

The underlying conditional distributions are characterized by the average {r)g i and variance af = VSf 



P(r\Si) = exp[-(r - {r)Y/2VS? }/^2nVSf . (Sll) 
which are each parameterized by the characteristic collaboration size Si- In cases where the average change (r) rj 0, 



then the distribution P(r\Si) is characterized by only the fluctuation scale di{r). Fig. S5 demonstrates that the 
normalized production change r-(i) = (r — (r;))/^ is distributed according to a Gaussian distribution. Hence, using 
normalized variables, we have mapped the process to a universal scaling distribution P(r\Si). 
When the distribution P(Si) is exponential, 

P(S l ) = Ae~ AS ' (S12) 
then mixture is termed an "exponential mixture of Gaussians" [44] . where the units have characteristic size Si — 



1/A. Fig. S10 shows that the distribution of collaboration radius Si is approximately exponential for each dataset, 
supporting the case for exponential mixing. Using the cumulative distribution of S for each data set we calculate 
A = 0.15 ±0.01 [A], A = 0.11±0.01 [B], and A = 0.11±0.01 [C]. While the tail behavior of P(r) can be used to better 
discriminate the value of ip, we do not have sufficient data in this analysis to perform a more rigorous test of the tail 
dependencies, or in general, to investigate the distribution of significantly large ri(t) values. 

The scaling relation <Xj(r) ~ Sf determines the functional form of the aggregate P^p{r). Clearly, er(r) increases 
for if) > values, whereas for values tp < 0, a(r) decreases with size S{. This latter case is empirically observed for 
countries and firms , whereby in general, large economic entities are able to decrease growth volatility by increasing 
and diversifying their portfolio of growth products. In our analysis of scientific careers we define Si = Med[ki(t)], the 
median number of distinct coauthors per year, as a proxy for the ability of the career to attract new opportunities, 
and hence, as a proxy for the size Si of an academic career. For professional athletes, we define the career size as the 
average number of points scored over the career Si = {Pi(t}). In Fig. H we calculate ip/2 w 0.40 ± 0.03 ( regression 
coefficient R = 0.77) for dataset [A], tp/2 w 0.22 ± 0.04 (R = 0.51) [B], and ip/2 « 0.26 ± 0.05 (R = 0.45) [C]. 

The role of mental, physical, and group spillovers is quite different in professional sports. Athletes attract future 
opportunities largely through their historical track record, which is heavily weighted on performance in the near past, 
and less on the cumulative history. Hence, for this performance-based labor force, we use a simple definition of "team 
value" to define the career size Si. This quantity is easier to define for basketball, since there are smaller differences 
between players of different team position than in other sports. For NBA player i we define Si as the average number of 
points scored per year, Si = (pi). Fig. S9 shows a crossover value S c which we interpret to reflect the fact that sports 
players typically fall into one of two categories: starters (everyday players) and replacement (game filler) players. We 
calculate ip/2 ss 0.38 ± 0.02 for emerging and "second string" careers with Si < S c , and a decreasing size variance 
relation (tjj < 0) for high- value careers with Si > S c . Similar values occur in the MLB. These two ip regimes reflect 
the crucial balance of risk and reward in short-term contract professions. 
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A variety of pdfs P^ (r) can result from the exponential mixture of Gaussians 

f°° 1 

P*(r)= Xe- xs —==eM-r 2 /2a 2 (r)}dS (S13) 
Jo y/2na z {r) 

depending on the value of if) which quantifies the size- variance relation. The functional form of P^{r) can vary in 
both the bulk and the tails of the distribution [33]. A simple result which follows from the case if) = 1 is the Laplace 
(double-exponential) distribution 



P^ =1 (r) = y — exp[-^ — Hj . (S14) 

This distribution is a member of the family of Exponential power distributions which follow from the range of values 
if> > [33]- In general, if the scaling values are in the range ip > 0, then the exponential mixture leads to an 
Exponential power distribution 

PW = ^w) exp[ -^ (w/ ' ),!| (S15) 

with shape parameter j3 in the range f3 € (0,2] [44]. The pure exponential P(r) with /3 = 1 corresponds to the case 
if) = 1. The pure Gaussian P(r) with /3 = 2 corresponds to the case if) = 0. 

Furthermore, if the annual production is logarithmically related to an underlying production potential, rii(t) oc 
lnUi(t), then r,(t) oc ln{/j(i) — laUi(t — 1) quantifies the logarithmic change ("growth rate") of Ui(t). This forms 
the analogy with growth dynamics of large institutions with size S ^> 1. For example, in the case of financial 
securities such as the stock of a company i, the growth rate measure the logarithmic change in the market's 
expectations of the company's future earnings potential captured by the market capitalization and price [SO]- As a 
result, distributions P(r) of career growth fluctuation r, which we plot in Figs.[2](a-d), can be seen as a bridge between 
the micro level and the macro level of economic growth fluctuation. A theory of micro growth processes can help 
improve the growth forecasts for economic organizations ranging in size from scientific collaborations to universities 
and firms [3J [3J 131 SZHSD]- 



IV. NONLINEAR PREFERENTIAL CAPTURE MODEL 



Here we describe a stochastic system in which a finite number of opportunities are distributed to a system of 
individual competing agents i = The opportunities are distributed in batches of P opportunities per arbitrary 

time interval. This model has two parameters. 

(i) 7r determines the preferential capture mechanism (the value it = 1 corresponds to the traditional "linear" 
preferential attachment model) and 

(ii) c determines the performance timescale 1 /c which is incorporated into the calculation of the capture rates of 
each individual. The value c = corresponds to a long-term memory and c ^> 1 corresponds to short-term memory. 

We use this simple model to show that a system governed by a preferential capture can become dominated by 
fluctuations when c is large. The value 1/c quantifies the "performance appraisal timescale": a small c corresponds 
to a labor system with long contracts, or some alternative mechanism that provides employment insurance through 
periods of low production, so that the ability to attract future opportunities is largely based on the cumulative record 
of career achievement. Conversely, a large c corresponds to a labor system with short contracts in which the ability to 
attract future opportunities is largely based on the accomplishments in the near past, requiring an agent to maintain 
relatively high levels of production in order to survive. In this latter case, we find that (natural) fluctuations in the 
annual production can cause a significant fraction of the careers to "fizzle out" leaving behind only a few "super 
careers" who attract almost all of the opportunities. In other words, short contracts can tip the level of competition 
into dangerous territory whereby careers are largely determined by fluctuations and not persistence. 



A. System of competing agents 

1) The system consists of I = 1000 agents competing for P opportunities that are allocated in a single period. 
There is no entry, hence the number / is kept constant. Also, P is also kept constant, so there is no growth in 
the labor supply. 
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2) We run the Monte Carlo (MC) simulation for T ee 100 time periods and all agents are by construction from the 
same age cohort (born at same time). 

3) Each time period corresponds to the allocation of P = E<=i "cm opportunities, sequentially one at a time, to 
randomly assigned agents i, where uq^ ee 1 is the potential production capacity of a given individual. 

4) The assignment of a given opportunity is proportional to the time-dependent weight (capture rate) Wi(t) of each 
agent. Hence, the assignment of 1 opportunity to agent i at period t results in the production (achievement) 
rii(t) to increase by one unit: rii(t) — > rii(t) + 1. In the next time period t + 1, we update the weight Wi(t + 1) 
to include the performance rij(i) in the current period. 



B. Initial Condition 



The initial weight at the beginning of the simulation is Wi(t = 0) = n c for each agent i with n c = 1. The value 
n c > ensures that competitors begin with a non-zero production potential, and corresponds to a homogenous system 
where all agents begin with the same production capacity. Hence, we do not analyze the more complicated model 
wherein external factors (i.e. collaboration factors) can result in a heterogeneous production capacity across scientists. 
By construction, each agent begins with one unit of achievement rii(t = 1) = 1. 



C. System Dynamics 

1) In each Monte Carlo step we allocate one opportunity to a randomly chosen individual i so that rii(t) — > rii(t) + l 

2) The individual i is chosen with probability Vi(t) proportional to [tu^t)]* 

vm = ™ i(t) ' (si6) 

where the value Wi(t) is given by an exponentially weighted sum over the entire achievement history 

t-i 

Wi(t) ee n *(* ~ Ai)e~ cAt . (S17) 

At=l 

The parameter c > is a memory parameter which determines how the record of accomplishments in the past 
affect the ability to obtain new opportunities in the current period, and therefore, the future. The limit c = 
rewards long-term accomplishment by equally weighting the entire history of accomplishments. Conversely, when 
c> 1 the value of Wi(t) is largely dominated by the performance n^t—1) in the previous period, corresponding 
to increased emphasis on short-term accomplishment in the immediate past. Intermediate values < c < 1 
weight more equally the immediate past and the entire history of accomplishment. 

3) The exponent n determines how the relative ability to attract opportunities Vi/Vj — [wi(t)/wj(t)] 7T depends 
on the weights Wi(t) and Wj(t) between two individuals i and j. The linear capture case follows from ir = 1, 
uniform capture it = 0, super linear capture 7r > 1, and sub- linear capture 7r < 1. 

4) At the end of each time period, the weight Wi (i) is recalculated and used for the entirety of the next MC time 
period corresponding to the allocation of the next I x n c achievement opportunities. 



D. Model Results 



We simulate this system for a realistic labor force size / = 1000 with the assumption that in any given period, 
an individual has the capacity for one unit of production (n c = 1). We evolve the system for T — 100 periods 
corresponding to / x n c x T Monte Carlo time steps. The timescale T represents the (production) lifetime of individuals 
with finite longevity. In this model we do not include exogenous shocks (career hazards) that can result in career 
death |16j . Here we analyze four quantities: 

1) The distribution P(N) of the total number of opportunities A^(T) ee Y^t=i n i(t) captured by agent i over the 
course of the T— period simulation. 



2) The distribution P(a) of the career trajectory scaling exponent on defined in Eq. S7 which quantifies the 
(de) acceleration of production over the course of the career. 

3) The distribution P{r) of production outcome change r defined in Eq. [2] which quantifies the size of endogenous 
production shocks. 

4) The distribution P(L) of career length Li which measures the active production period of each career starting 
from t — 0. We define activity as the largest period value Li for which rii(Li) = 0, which in other words, 
corresponds to truncating all production values from the end of the trajectory rii(t) and defining Li as the 
length of this time series. 

We display these four distributions, from left to right, for varying it and c values, in each panel of Figs. |S12| 



S16 Empirical distributions calculated from MC simulations are plotted as blue dots, with benchmark distributions 
described below plotted as solid green curves. For each tt and c value we simulate 10 MC systems, and combine the 
results into aggregate distributions which are shown. For simulations with tt > 1 the pdf data are aggregated over 
the results of 50 MC simulations. We list below some of our main observations. 

For 7r = 1, independent of c, we observe exponential P(N), consistent with the prediction of the linear preferential 
capture model in the case of no firm entry (b = 0) in the model of Kazuko ct al. [43 . However, the distribution P(L) 
and the distribution P(a) does depend strongly on c, reflecting the possibility of career "sudden death" for large c. 

For the P(a) distributions (middle-left panels), the solid green line is a best-fit Gaussian distribution (using the 
MLE method) for the set of values computed for careers that did not undergo "sudden death." 

For the P(r) distributions (middle- right panels), the solid green curve corresponds to a best-fit Laplace distribution 
(using the MLE method) and the dashed red curve corresponds to a best-fit Guassian distribution (using the MLE 
method) which we show only for benchmark comparison. Typical empirical distributions (values shown as blue dots) 
range from being distributions that are Gaussian to distributions that are Laplacian in the bulk but with heavy tails. 

For the P(L) distributions (right most panels), we note that the most likely career length L is typically either 
L = 1 or L = T for all systems analyzed. However, there are likely c and tt parameter values corresponding to P(L) 
that is uniform distributed over the entire range of L values, which may be an interesting class of system to analyze 
in future analyses since such a system promotes diversity across the entire longevity spectrum. The system we show 
for 7r = 1.2 and c = 1 appears to be close to this scenario. 



Fig. |S12| shows the null model with no preferential capture (it — 0). We confirm that the careers in this model 
are driven by a stochastic accumulation process that is equivalent to a Poisson process with rate A p = 1. In this 
homogenous system, each career gains on average one opportunity each time period, so that at the end of the 
simulation, the distribution P(N) is a Poisson distribution with (N) = X p T (shown as the solid blue line) which fits 
the model data excellently. For these careers, the typical a = 1, the production changes are well-approximated by 
a Gaussian distribution, and most careers are sustained for the maximum possible lifetime corresponding to T periods. 



Fig. S13 shows the system with c = corresponding to comprehensive career appraisal corresponding to a 
long-term memory system. We analyze this system for 4 values of it = 0.8,1.0,1.2,1.4. This "long-term memory" 
scenario corresponds to a long-term contract profession whereby careers are less vulnerable to periods of low 
production. As a result, most careers sustain production throughout the career. 



Fig. S14 shows the system with c = 0.1 corresponding to an effective memory timescale of 1/c = 10 periods. We 
analyze this system for 4 values of it — 0.8, 1.0, 1.2, 1.4. This "medium-term memory" scenario yields a rich variety 
of careers for it = 1, but for n = 1.2 the system becomes quickly dominated by "rich-get-richer" effects which results 
in careers being vulnerable to low production fluctuations. 



Fig. S15 shows the system with c = 1 corresponding to an effective memory timescale of 1/c = 1 period. We 
analyze this system for 4 values of it = 0.8,0.9,1.0,1.1. For all values of it analyzed, we observe a system that is 
dominated by careers that are cut short by the high levels of competition induced by the relatively high value placed 
on continued production. 



Fig. S16 shows the extreme case of a "no memory" scenario in which i0j(t) ~ nj(f — 1) whereby most careers 



experience sudden death due to endogenous negative production shocks early in their career. The lucky few careers 
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who survive this period end up as rich-get-richer "superstars." This behavior occurs for all systems analyzed using 4 
values of tt = 0.8, 0.9, 1.0, 1.05. 



E. Discussion of the model in relation to the Academic labor market 

One serious drawback of short-term contracts are the tedious employment searches, which displace career momentum 
by taking focus energy away from the laboratory, diminishing the quality of administrative performance within 
the institution, and limiting the individual's time to serve the community through external outreach [2 [6]. These 
momentum displacements can directly transform into negative productivity shocks to scientific output. As a result, 
there may be increased pressure for individuals in short-term contracts to produce quantity over quality, which 
encourages the presentation of incomplete analysis and diminishes the incentives to perform sound science. These 
changing features may precipitate in a "tragedy of the scientific commons." 

Aside from promoting circumspect research, job security in academia diminishes the incentives for scientists to 
"save and store" their knowledge for future liquidation in the case of employment emergency, and thus promotes the 
institution of "open science" PQ. However, a policy shift towards short-term contracts, along with the heightened value 
of intellectual property, may alter the course of publicly funded "open science." This scientific commons emerged 
from the noble courts during the Renaissance as a hallmark of the scientific revolution and now faces pressure from 
what has been termed "intellectual capitalism," with the vast privatization of knowledge and innovation ("closed 
science") occurring in public universities and corporate R&D lj. An academic system that is dominated by short 
term contracts, stymied by production incentives that favor quantity over quality, and jeopardized at the level of the 
"open knowledge" commons, presents a new institutional scenario revealing selection pressures that could alter the 
birth and death rates of high-impact careers. 

The purpose of this stochastic model is to show how careers can become very susceptible to negative production 
shocks if the labor market is driven by a preferential capture mechanism with 7 > 1 whereby early success of an 
individual can lead to future advantage. However, this model also shows that the onset of a fluctuation-dominant 
(volatile) labor market can also be amplified when the labor market is governed by short-term contracts reinforced by a 
short-term appraisal system. In such a system, career sustainability relics on continued recent short-term production, 
which can encourage rapid publication of low-quality science. In professions where there is a high level of competition 
for employment, bottlenecks form whereby most careers stagnate and fail to rise above an initial achievement barrier. 
Instead, these careers stagnate, and in a profession that shows no mercy for production lulls, these careers undergo a 
"sudden death" because they were "frozen out" by a labor market that did not provide insurance against endogenous 
fluctuations. Such a system is an employment "death trap" whereby most careers stagnate and "flat-line" at zero 
production. However, at the same time, a small fraction of the population overcomes the initial selection barrier and 
are championed as the "big winners" , possibly only due to random chance. 

Table demonstrates how the life expectancy decreases with increasing c even for the linear preferential capture model 
corresponding to tt — 1. With increasing c, the model simulates systems with shorter contracts (shorter appraisal 
"memory" timescales), and so larger percentages of the population die before characteristic ages T c (p), values that 
decrease with increasing c for a given p. 
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T c (p) as a % of T, (% T) 



p = 0.1 


p = 0.25 


p = 0.5 


p = 0.75 


c 


U (long term) 


0.94T 


0.98T 


LOOT 


LOOT 


r 


0.1 


0.20T 


0.79T 


0.99T 


LOOT 


r 


1.0 


0.01T 


0.02T 


0.05T 


0.15T 


f 


10.0 (short term) 


0.01T 


0.01T 


0.02T 


0.06T 



TABLE SI: Decrease in career life expectancy as a result of short-term contract length in the tt = 1 linear preferential capture 
model. The fraction p of the population that experienced career termination before the crossover age T c (p): "p percent of the 
population died before reaching the age L — T c (p)." As c increases (recall the appraisal "memory" timescale is 1/c) towards a 
short-term contract scenario, a significant fraction of the population (increasing p) dies before reaching a smaller and smaller 
T c (p). The empirical value of T c (p) is given as a percentage of the maximum career length T corresponding to the stopping time 
of the Monte Carlo simulation. The value T c (p) is calculated using the equality p — CDF(T < T c (p)), where CDF(T < L) 
is the cumulative distribution function of career length L. To estimate CDF(T < L), we combine an ensemble of 10 MC 
simulations for each c value. In the model simulations we use T = 100 periods. 




FIG. SI: Positive career shocks likely associated with reputation boosts. Examples of career production trajectories Ni(t) 
that have significant deviations from the scaling hypothesis in Eq. |S7 | These significant deviations likely follow extraordinary 
scientific discoveries (and the publicity and reputation that are typically rewarded) which can vault a career and result in 
lasting benefits to the individual. 
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FIG. S2: Regularities in the career trajectory Ni(i). We analyze the normalized career trajectory iVj(t) = Ni(t)/{n,i) which 
allows us to aggregate Ni(t) with varying publication rates (rii). As a result, we can better quantify the scaling exponent 
a which quantifies the acceleration of the typical career over time. We calculate a using OLS regression on log-log scale of 
the average normalized career trajectory (N'(t)} = ^^y^- For reference, each AT|(i) trajectory in panels A, B, and C has a 

corresponding best-fit curve that is a dashed line. (A) For the scientific careers, we calculate a values: 1.28 ± 0.01 for Dataset 
A, 1.31±0.01 for Dataset B, and 1.15±0.02 for Dataset C. These values are all significantly greater than unity, 57 > 1, indicative 
of a systematic cumulative advantage effect in science. (B) The standard deviation aN'(t) has a broad peak, likely related to 
career shocks that can significantly alter the career trajectory. (C) The average normalized career trajectory for NBA careers 
has a ~ 1 (D) The average normalized career trajectory for MLB careers has a ~ 1. For visual comparison, the solid straight 
black line in panels A,B and C correspond to a linear function with a = 1. 
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FIG. S3: Using scaling methods to show approximate data collapse of each Ni(t). Normalized trajectory Ni(t) = Ni(t)/Ai 
plotted using the scaled time t' = t ai for each career over the time horizon t £ [1,40] years. We plot the 100 Ni(t) curves 
belonging to datasets [A], [B], and [C] in the corresponding panels. There is approximate data collapse of all the normalized 
trajectories Ni(t) along the dashed green line corresponding to the rescaled career trajectory Ni(i) = t' with a' = 1 by 
construction. We also plot in red the corresponding average value (Ni(t)) with la error bars for logarithmically spaced t' 
intervals. Deviations from (Ni(t)) are indicative of career shocks which can significantly alter the career trajectory. 



1.5 



0.5- 



2 - □ Dataset A 






0.8 1.2 1.6 2 

a. 



0.8 1.2 1.6 2 

a. 



0.8 1.2 1.6 2 

a. 



FIG. S4: Increasing returns to scale a > 1. Probability distribution of the individual Oi values calculated for each career 
using the scaling model 2Vj(i) ~ t ai over time horizon t £ [1,40] years. The average (o^) and standard deviation a(ai) for each 
dataset are: 1.42 ± 0.29 [A], 1.44 ± 0.26 [B], 1.30 ± 0.31 [C]. The distribution of oti values indicate that career trajectories are 
typically accelerating (a« > 1), most likely the result of a cumulative advantage effect. 
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FIG. S5: Universal patterns in underlying production fluctuations of scientists. Accounting for variable individual publication 
factors, such as academic subfield or group collaboration size, we find that the normalized annual production change r^(t) = 
[ri(t) — ( r )i\l<7i is distributed according to a Gaussian distribution, with (r'} = and cr(r') = 1 by construction (solid lines 
show best-fit Guassian distributions using the maximum likelihood estimator method). This results indicates that the Laplace 
distribution shown in Fig. [2]results from a mixture of Gaussian distributions Pi{r = air') that indicate that annual production 
is consistent with a proportional growth model.. 
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FIG. S6: Universal patterns in the production fluctuations of athletes. For athlete careers in the NBA and MLB we define 
production change for (A,C) the change in the number of in-game opportunities and (B,D) the change in the number of in- 
game successes. (A,B) Since the detrended production change R is defined to have standard deviation a = 1, the pdfs P(R) 
approximately collapse onto a universal "tent-shaped" Laplace pdf (solid green line). (C,D) For sports careers, we also define 
a measure R' which account for variable individual production factors, such as propensity for injury, team position, etc. As a 
result normalized annual growth rate R'i = [zi(t) — (z(t))]/cz(t) is normalized twice, once to account for age factors and once 
to account for individual factors. The quantity Zi(t) = (ri(t) — (fi))/cr, is normalized with respect to individual factors, where 
(ri) and Oi are the average and standard deviation of the production change of career i. Then, we aggregate all Zi(t) values for 
a given career year t in order to calculate the average (z(t)) and standard deviation cr z (t) over an careers. The final quantity 
R'i represents a normalized annual production change which is distributed in the bulk according to a Gaussian distribution, 
with (R') « and cr(r') w 1 by construction (solid lines show best-fit Guassian distributions using the maximum likelihood 
estimator method). This results indicates that the tent-shaped distributions in (A,B) results from a mixture of conditional 
Gaussian distributions Pi(R = aR') that indicate that annual production is consistent with a proportional growth model. 
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FIG. S7: Universal micro-scale output distribution P(Q) which accounts for coauthorship variability. The normalized output 
Q oc rii/kj 1 is a residual output after we quantitatively account for the collaboration size ki corresponding to the number of 
distinct coauthors of author i. Each pdf is well-approximated by the Gamma distribution P(Q) oc Q m_1 exp[— Q/0] which 
suggests that production at the micro scale is governed by a Gamma Levy process. We calculate the Gamma distribution 
parameters using the maximum likelihood estimator method (distributions shown by solid and dashed curves), and find an 
insignificant difference between [A] and [B] scientists with Gamma shape parameter m and scale parameter 8. However, for 
dataset [C] scientists, the output distribution is more skewed towards smaller Q values, possibly reflecting the relative advantage 
that senior scientists gain due to reputation, experience, and knowledge spillover factors. 




FIG. S8: Aggregate production distributions can be deceiving. Unconditional distribution of annual publication rate n(t) 
appears as log-normal distributions because it is a mixture of underlying distributions that depend strongly on collaboration 
factors. We define rii{i) as the number of papers published in (A) At = 1 and (B) At = 2 year periods, which reduces the 
finite-size effects arising from the calendar year labeling of publication dates. (A) We combine m(t) values for all values of 
t, and find excellent agreement between the empirical P(n(t)) data points and the log-normal model. We use the maximum 
likelihood estimator method to calculate the log-normal parameters gl = a(lnn) and fj, = (Inn). (B) In order to analyze the 
time-dependence of P(n(t)), we separate rii(t) values from Dataset A into 5 subsets, depending on the range t years into the 
career, as indicated in the figure legend. We offset each pdf by a constant factor in order to distinguish each pdf, which are 
also well-approximated by log- normal distributions (shown as solid curves). 
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Average points scored per year, S. = <p > Average strikeouts per year, S. = < K. > 




Average hits per year, S. = < H > 



FIG. S9: Quantifying the growth fluctuations of sports careers. The size variance relation for sports careers is similar to 
academic careers for small Si. However, for relatively large Si the relation becomes decreasing corresponding to ip < 0, 
analogous to what is found for firm growth 1341 148H50] . The decreasing relation for Si > S c likely follows from the fact that 
in sports, there is a hard upper limit to the number of opportunities available to a player in a given year. Hence, individuals 
with large Si are likely the starters on their teams, since it is neither economical nor in the strategy of winning to keep players 
above a threshold value S c out of the game, and so these players typically remain as positional starters except for episodic 
leaves of absence due to injury. Hence, these players experience smaller <Ji(r) due to limitations to their potential for further 
career growth. However, players with Si < S c are typically on the fringe of being released or provide alternative value to the 
team, and so these individuals experience larger fluctuations in team play because they are easily dispensable, especially in a 
profession dominated by short-contracts lasting sometimes less than a year. For each dataset, we use careers with career length 
Li > 3 seasons. (A) NBA basketball players: Units of Ui{R) are normalized minutes played. We define the scaling relation 
ai(R) ~ <p l ) v ' /2 between the average number of points scored per season (jpi) = X!t=i Pi{t)/Li and the standard deviation (Ji(R). 
In this way, we utilize the average points per season as the proxy for the ability of a player to obtain future opportunities 
which are realized as minutes played. Using S c = 720 points, we calculate tp/2 = 0.38 ± 0.02 (regression coefficient R = 0.50 
and AN OVA F-test significance level p » 0) for Si < S c and tp/2 = -0.25 ± 0.07 (R = 0.15 and p w 10" 3 ) for S t > S c . 

(B) MLB pitchers: Units of Oi(R) are normalized IPO (innings pitched in outs). Interestingly, <7i(R) continues to increase 
for Si > S c , possibly due to the relatively high career risk attributed to throwing arm injury. Using S c = 65 strikeouts, we 
calculate tp/2 = 0.37 ± 0.01 (R = 0.48 and p « 0) for Si < S c and tp/2 = +0.15 ± 0.07 (R = 0.07 and p w 0.02) for Si > S c . 

(C) MLB batters: Units of a t {R) are normalized AB (at bats). Using S c = 68 hits, we calculate ip/2 = 0.44 ± 0.01 (R = 0.59 
and pRiO) for S t < S c and tp/2 = -0.37 ± 0.03 (R = 0.21 and p w 0) for Si > S c . The dashed black (blue) line in each panel 
is a least squares linear regression on log- log scale for all data values with Si less (greater) than S c . The data shown with error 
bars represent the average (<7j(i?)) and corresponding 1 standard deviation values calculated using equally spaced Si bins on 
the logarithmic scale. 
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FIG. S10: Exponential distributions of coauthor radius in Physics. We test the hypothesis that the distributions P(r) 
for annual production change r (shown in Fig. [2| follow from an exponential mixing of Gaussians with varying fluctuation 
scale Oi oc Med[ki(t)]^^ 2 . An important criteria for this model is that the distribution of Si = Med[ki(t)\ is exponential, 
P{Si) ~ exp[— XSi\. We plot the cumulative distribution function (CDF) P(x > Si) for each dataset, and confirm that the 
distributions are approximately linear on log-linear axes. Using linear regression, we calculate A = 0.15±0.01 [A], A = 0.11±0.01 
[B], and A = 0.11 ±0.01 [C]. 
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FIG. Sll: Approximately exponential distribution of scoring value in the NBA. We further test the hypothesis that the 
distributions P(R) for annual production change R in professional sports (shown in Fig. [2]C and D) follow from an exponential 
mixing of Gaussians with varying fluctuation scale ai oc (pi)^^ 2 . An important criteria for this model is that the distribution 
of "team value" is exponential, P({Pi)) ~ exp[— A(pj)]. We plot the cumulative distribution function (CDF) P(x > {pi}) 
for each dataset, and confirm that the distributions are approximately linear on log-linear axes. We show the CDFs calculated 
using all careers with career length Li > L c years, for L c = 1, 3 years. 
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FIG. S12: A production output null model with tv = agrees with the predictions of a Poisson process. (Far left) The 
cumulative distribution C'DF(x > TV) is in excellent agreement with the prediction of a Poisson process with rate X p — 1 and 
corresponding average (N) = X P T = 100. The solid green curve is the corresponding Poisson CDF using (TV) = 100. (Middle 
left) Furthermore, the typical scaling exponent (a) — 1 which is also consistent with Poisson trajectories. (Middle right) The 
distribution of production changes is close to Gaussian. (Far right) The typical career length Li spans the entire system length 
T, indicating low levels of career risk. 
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FIG. S13: The production output model with c = 0. Results of MC simulations for a "long-term appraisal" scenario. Careers 
are less vulnerable to low-production phases, and as a result, most agents sustain production throughout the career for a 
relatively large range of it values. 
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FIG. S14: The production output model with c = 0.1. Results of MC simulations for a "medium-term appraisal" scenario. 
The corresponding memory time scale is approximately 10 time periods, and so only for significantly large ?r = 1.4 do we 
observe a labor market scenario in which there is a significant death rate and just a few "big winners" corresponding to those 
agents with a > 1. 
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FIG. S15: The production output model with c = 1.0. Results of MC simulations for a "short-term appraisal" scenario. The 
corresponding memory time scale is approximately 1 time period. Even for ir < 1, the system is driven by fluctuations that can 
cause career "sudden death" for a large fraction of the population. For n > 1 we observe a very quick transition to a significant 
death rate and just a few "big winners" corresponding to those agents with a > 1. 
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FIG. S16: The production output model with c = 10.0. Results of MC simulations for a "zero-memory appraisal" scenario 
wherein only the previous period matters, Wi(t) = rii(t — 1). Even for linear preferential capture 7T = 1, the systems shows "no 
mercy" for careers that are stagnant for possibly just one period. As a result, just a few "lucky" agents are able to survive 
the initial fluctuations and end up dominating the system. For n values close to unity, n — >■ 1, the systems quickly becomes an 
employment "death trap" whereby most careers stagnate and "flat-line." 



